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COMPUTE  GENERATION  OF  RANDOM  NUMBERS 


x  ,  INTRODUCTION 


A  random  number  is  a  quantity  which  is  specified  probabilistically, 
and  a  random  number  generator  is  a  device  which  outputs  numbers  under  the 
constraints  of  probability  lavs-  The  probability  laws  governing  the 
sequence  describe  (l)  the  statistical  properties  of  the  output  sample, 
and  (2)  the  relationship  between  successive  numbers  in  the  sequence.  The 
present  memorandum  reports  progress  in  the  development  of  computer  tech¬ 
niques  for  generating  randan  numbers. 


The  probability  law  governing  the  output  sample  often  takes  the  form 
of  a  theoretical  probability  distribution  function  which  the  output  as  a 
whole  is  expected  to  satisfy.  The  usual  sequential  requirement  is  that 
successive  numbers  in  the  output  sequence  be  statistically  independent. 

There  are,  of  course,  numerous  possible  rules  to  govern  non-independent 
sequences. 

V-x 

Random  numbers  have  their  greatest  use  in  simulation.  A  random 
number  generator  is,  in  fact,  a  program  to  simulate  sample  distributions. 

A  simulated  sample  is  often  very  useful  in  solving  psychometric  problems, 
because  its  properties  can  be  controlled  more  easily  than  can  those  of 
empirical  samples.  (See,  for  example,  the  study  of  an  allocation  tech¬ 
nique  applied  to  current  aptitude  input  by  Boldt,  Wiskoff,  and  Fitcb,  i960. ) 
However,  the  real  utility  of  random  number  generators  is  achieved  in  simu¬ 
lation  of  theoretical  systems  and  probabilistic  control  processes.  A 
typical  theoretical  system  often  studied  by  simulation  is  a  mathematical 
model  of  behavior.  Reasonable  models  of  behavior  often  turn  out  to  be 
mathematically  intractable.  However,  canputer  programs  may  be  written 
which  simulate  the  functional  operations  of  the  model.  It  would  then  be 
possible,  with  a  Judicious  choice  of  parameters,  to  run  many  empirical 
simulations  of  the  mode 1 .  In  this  way,  the  model  may  be  tested  so 
thoroughly  that  explicit  mathematical  solutions  of  the  model  equations 
are  no  longer  needed.  Random  number  generators  are  essential  in  this 
endeavor,  because  most  psychological  models  are  written  in  the  language 
of  probabilities. 


An  example  of  a  probabilistic  control  process  is  a  psychological  ex¬ 
periment  in  which  the  stimulus  presentations  (or  the  administered  rewards) 
have  sane  chance  dependence  on  other  events.  Experiments  sometimes  re¬ 
quire  that  stimuli  be  chosen  "at  randan",  or  with  some  probabilistic 
dependence  on  the  subject's  last  response,  or  at  certain  "randan"  intervals 
in  time.  For  example,  studies  of  multi-person  interaction  often  use  a 
complicated  reward  structure  in  which  the  outcome  administered  to  any  in¬ 
dividual  depends  on  the  Joint  response  pattern  of  the  group  as  well  as 
upon  chance  events  set  up  by  the  experimenter.  The  probabilistic  character 
of  these  contingencies  can  often  be  achieved  by  a  randan  number  program. 

In  fact,  there  are  many  interesting  experiments  which  cannot  be  done  with¬ 
out  the  canputer' s  probabilistic  control. 


METHODS 


The  core  of  most  random  number  programs  is  a  simple ,  repeat-able 
arithmetic  operation.  In  almost  every  program,  this  core  operation  pro¬ 
duces  numbers  which  are  uniformly  distributed  along  tne  interval  from  zero 
to  one.  This  is  by  no  means  a  necessary  first  step,  but  it  makes  possible 
a  variety  of  second  steps,  such  as  converting  uniformly  distributed  numbers 
into  other  distributions.  For  example,  it  is  much  easier  to  convert 
uniform  random  variables  into,  say,  a  beta  distribution  than  to  construct 
the  beta  distribution  directly.  Therefore,  most  random  number  programs 
will  have  two  parts :  a  subroutine  to  generate  uniform  random  numbers  and 
a  master  program  to  make  them  conform  to  same  probability  distribution. 

Of  the  many  computer  metnods  for  generating  uniformly  distributed 
numbers,  none  are  perfect.  All  produce  a  repeating  finite  sequence,  the 
terms  of  which  can  be  perfs  bly  predicted  from,  '-ha  initial  conditions. 

Nevertheless,  the  sequences  may  be  extremely  long  (so  that,  for  practical 
purposes,  they  do  not  repeat),  and  they  may  possess  most  of  the  statistical 
properties  of  a  truly  random  sequence.  For  this  reason,  the  numbers  pro¬ 
duced  are  sometimes  called  "pseudo-random".  A  discussion  of  alternative 
methods  for  generating  uniform  random  numbers  is  beyond  t.he  scope  of  the 
present  memorandum,  but  the  one  most  canmonly  used  will  be  briefly 
described. ^ 


THE  POWER  RESIDUE  METHOD 

An  arithmetic  process  that  produces  a  long,  erratic  sequence  can  be 
invented  easily;  to  make  the  terms  independent  and  equi-probable  is  much 
more  difficult.  The  power  residue  method  canes  close  to  satisfying  all 
these  requirements.  It  relies  on  computations  of  the  residues  (remainders) 
of  successive  powers  of  a  number,  x  modulcr^M,  where  M  is  large  compared 
to  x.  In  particular,  if  x  and  M  are  relatively  prime^S  the  sequence 

i  =  x*1  (mod  M),  or  the  sequence  u  =  a.xr‘  (mod  M)  where  a  is  also  rela- 

.  *4 

tlvely  prime  to  M,  is  a  repeating  sequence  with  random  properties. 


^  Alternate  methods  are  given  in  Hamming  (1962);  IBM  Reference  Manual 
C20-8011;  and  Meyer  (1956). 

^The  expression  u  =  v  (mod  M)  means  that  (u  -  v)  is  divisible  by  M.  The 
operation  u.v  (mod  M)  is  read  "multiply  u  by  v,  divide  the  product  by  M, 
and  take  the  least  positive  remainder". 

^Two  numbers,  u  and  v,  are  relatively  prime  if  their  greatest  common 
divisor  is  one. 


The  residue  sequence  depends  on  the  relationship  between  x  and  M.  if 

M  =  2^  for  same  positive  int-ger  b,  and  if  x  is  a  suitably  chosen  odd 

integer,  there  will  be  2b  *  terms  produced  before  t.ne  sequence  repeats. 

This  fact  oas  s  me  practical  Importance.  To  make  the  sequence  as  long  as 

possible  in  a  bleary  compu*  Log  machine,  it  Is  convenient  to  choose  M  to 

be  the  maximal  word.  If  there  are  b  "tits"  available,  the  maximal 
b  19 

computer  word  is  2  (°r.  the  GE  225,  M  =  2  ,  and  the  maximal  sequence 

17  • 

length  is  t  ere  fore  2  .  >  Tils  choice  ras  the  virtue  +  hst  the  dlvlsloh 


operation  ir  ue  formaVon  of  residues  is  eliminated 
multiplies  xr‘  by  xr‘  ,  tie  product  "overflows"  (i.e 


when  the  c computer 
,  xn  +  (n  ‘  1}  >  219), 


and  is  aut ■.•man  a^iy  ^<-r  r  ,uau  :  tv?  resivue. 


Tae  power  residue  awr*.‘::d  nas  t:o  programmed,  tested,  and  ir.corpo 
rated  in  iara-  r&nd.m  um*.c-r  programs.  Three  versions  exist,  one  written 
in  GE  225  ms  -Lae  ^ar.guag-  'jAtj,  one  written  in  FORTRAN  for  the  GE  225, 
and  oh“  wr itf®"  la  na  tli-  ja.r.guage  for  the  IBM  7090  A  statistical 
analysis  of  *  :.e  resul"?  ha?  been  completed  and  will  be  reported  separately. 
All  toe  pr  grams  perform  wen  or  tests  of  uniformity  and  independence. 


METHODo  FOR  CONVERSION 

Ltiept  r.de  <■  *  ur  If  nrii.v  d  Ler-ituted  random  numbers  are  rarely  all  that 
is  needed  ’ c  scive  a  pr'^lem.  viost  simuiati'ons  require  that  the  numbers 
sa'isfy  some  additional  trofcab ility  lav.  Therefore,  general  methods  are 
needed  for  conve  f  :g  uniform  random  numbers  into  other  more  complex  prob - 
r.blii*  y  distributions.  Two  conversion  techniques  can  be  distinguished: 

(i)  integration  neVio&e  and  (.-  )  composition  methods. 

Integration  Methods.  Let  f(x)  be  a  continuous  probability  density 
function  which  is  being  simulated  and  let  {u)  be  a  sequence  of  uniformly 
distributed  random.  numbers  on  the  unit  interval.  The  u's  are  interpreted 
e*  percentile  scores.  For  each  value  of  u,  we  wish  to  find  a  number,  v, 
which  is  at  the  uth  percentile  of  the  distribution  f(x).  The  v's  will 
then  have  the  required  f(x)  distribution.  For  example,  if  f(x)  is  tie 
standard  normal  distribution,  and  if  u  is  0.975,  then  the  value  of  v  which 

is  greater  than  exactly  97^  of  the  normal  distribution  turns  out  to  be 

v  =  1  96.  The  value  1.96  can  be  found  in  the  usual  tables  of  the  normal 
distribution  (although,  of  course,  tables  would  not  be  used  with  the  random 
number  generator).  Mathematically,  The  example  from  tne  normal  distribu¬ 
tion  is  expressed  as  follows:  If 

v  v  -1 

[’  ll  -x*/2 

u  *  0.975  =  f(x)dx  =  ^2n)  e  dx , 
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then  v  =  I.96.  More  generally,  given  u  one  mush  be  able  to  find  v  such 
that  the  area  under  f'„x. fr:*n  ite  lever  limit  to  v,  equals  u: 

V 

r 

u  =  J  f  ’x  ''ix  [l] 

—CO 

It  is  this  integration  problem  -.rat  gives  the  method  its  name. 

Whether  v  is  easily  determined  depends  on  the  density  function  f . 
Soros  functions,  such  as  the  exponential, 

f;x)  =  Xe"**  ,  X  >  0  [2] 


or  the  Laplace, 

f(x)  *  ~  Xe  ,  X  >  0  [j] 

are  easily  integrated,  and  an  exact  solution  to  equation  [l]  can  be  pro 
grammed.  This  programming  has  beer,  done  for  the  exponential  and  the 
Laplace  distributions.  Samples  from  these  random  number  routines  can  be 
quickly  generated.  Because  their  accuracy  depends  only  on  the  residue 
calculations,  the  samples  conform  well  to  taeoretical  distributions.  The 
exponential  distribution  occurs  often  in  models  of  latency  mechanisms  (e.g. , 
models  for  visual  search  time  or  detection  response  time).  It  is  also  a 
component  of  more  complex  distributions,  such  as  the  gamma,  and  is  likely 
to  be  one  of  the  most  useful  random  number  programs. 

The  generality  of  the  integration  method  depends  on  hew  easily  Equa¬ 
tion  [l]  can  be  solved.  When  ar.  exact  solution  cannot,  be  found,  same 
approx ima-cion  technique  is  often  useful  This  is  the  case  when  f  is  a 
normal  distribution.  The  most  convenient  method  for  integrating  the  normal 
distribution  on  a  computer  is  to  use  an  approximation  formula  (Hamming, 

1962;  Ralston  and  Wilf,  1962).  Usually,  use  of  this  formula  entails 
calculating  the  first,  few  terms  of  an  infinite  series.  The  normal  distri¬ 
bution  subroutine — a  FORTRAN  pr ogram  for  the  GE  225  or  IBM  7090 — calculates 
the  area  under  any  continuous  segment  of  the  normal  curve,  using  an 
approximation  formula.  This  program  is  quite  fast,  as  it  must  he  when 
large  samples  are  needed,  and  it  is  accurate :  the  output  agrees  with  the 
tabled  normal  distribution  to  six  decimal  places. 
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When  nel'  *-,er  ar.  acaly':  soiu*  lot 


available  for  Equation  [’]  *  '  *  '“wa’ 

iterative  procedure,  Several  iterative 
converting  uniform  random  rs  .  ;ic  a 


•  an  appro*  cna* e  formula  is 
rfly  b®  '•arr'i®d  out  by  an 
•e  nlquee  w -  ■  f-  investigated 
gamma  dietr ibur ion. 


for 


f .  x 


X  „  r  -1  -Xx 
=•■  c"  Xx  e 
1  r  • 


[4] 


where  X  >  0,  Tir)  =  J0  x1'1"  e"x  dx .  and  r  >  0.  All  ♦ae  iterative  met  beds 

use  a  series  of  quadratures  of  the  area  under  the  probability  curve.  At 
each  step,  an  ordinate  grid  is  overlayei  on  'he  area  'o  be  Integrated. 

This  grid  divides  the  area  into  four-sided  geometric  figures  defined  by 
two  adjaee- 1  ordinals  or  th*-  str'd  *he  «ea»®n+  of  the  x  axis  between  the 
ordinat.es,  and  a  line  segment  near  t  :e  probability  curve  chosen  to  simplify 
computations .  Tap  total  ar~a  under  -"urve  is  approximated  by  the  sum 
of  the  small  grid  areas.  Figure  1  illustra'es  this  method.  On  earn 
iteration,  the  number  of  grid  lines  is  Increased,  making  the  approx imat Ion 
closer  to  the  true  area  until  two  successive  iterations  give  estimates 
which  differ  by  less  than  a  pre  set  tolerance.  In  tie  FORTRAN  program 
that  uses  the  method  shew,  in  Figure  1,  N  rectangles  are  used  to  determine 


N- 1 


s-Sfl  f  [K^r)] 


f(x)  dx 


t5] 


i^O 


k  -1 

on  the  first  iteration,  and  N*2  are  used  on  the  kth  iteration. 

The  most  time-consuming  step  in  computing  Equation  [5]  is  cal  cilia' ion 
of  the  ordinates,  f(i).  For  the  gamma  distribution,  Equation  ^ ,  the  GE  225 
computer  used  about  1*5  minutes  to  compute  1000  values  of  f(t).  The  number 
required  to  reach  a  tolerance  of  1 $  varied  from  about  100  ordinates  to 
10,000  ordinates,  depending  on  the  parameter  r.  Iterative  routines  that 
used  Simpson's  rule  or  trapezpoidal  approximations  were  equally  time  con¬ 
suming.^  These  results  force  the  conclusion  that  although  iterative 
routines  are  useful  for  integration,  they  are  not  fast  enough  to  provide 
random  numbers  in  quantity.  Methods  are  needed  which  will  not-  require  as 
many  operations.  Monte  Carlo  methods  may  satisfy  this  requirement. 

A  Monte  Carlo  method  is  a  procedure  by  which  random  processes  simulate 
or  approximate  mathematical  functions.  The  best-known  use  of  random  proc¬ 
esses  in  computers  is  Monte  Carlo  integration:,  to  find  the  area  under  a 


^Ralston  and  Wilf  (1^62)  give  more  thorough  presentations  of  iterative 
quadrature  methods. 
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carve,  the  computer  sca**.-'s  pc  i**.=  a-.  rard-m  in  a  space  of  known  area 
which  includes  tie  reglo:  to  h  e  integrated  as  a  props*  *  it  space.  The 
proportion  of  points  falling  :  side  '  ir *e fixation  rsa:or  is  the  estimate 
of  Its  area  To  genera*?  prcfcablii*  /  disc  ..hut ions  we  use  a  variation 
of  this  technique  attributed  to  I.  V  •  Neumann.  Pairs  of  independent 
uniform  random  numbers  (x,y)  are  ger.-rar.ed  x  In  tie  data  la  of  tie  prob¬ 
ability  density  function  f,  at  i  y  in  Its  range.  Each  x  a  "candidate" 
for  inclusion  in  tie  sane  1  ?  of  random  mashers  having  t.c.e  f  x)  distribu¬ 
tion.  If,  for  a  given  x,  f  ,x )  >  v  x  is  "accepted"  in  t..e  sample, 
otherwise  x  is  discarded,  ana  a  *ew  \x  v  pair  is  generated.  Figur-  2 
lllustra*es  this  procedure.  Ice  proport' r.*-  of  x  values  included  in  the 
sample  depends  on  tie  state  Qf  t  If  *  -  function  is  tlgily  peaked,  its 
range  will  he  large  compared  *■.  -  prcbatl2-.lt/  ordinates  at  most  values 

of  x,  and.  because  a  large  percentage  of  y  values  will  Men  exceed  the 
ordina+es,  many  x  value?  will  he  it- carded »  The  efficiency  of  the  me* rod 
depends  or  1 1-  probatil'-  ‘  .a*  x  •  s.  o:  discarded,  w.-.  l cl  is 


p  { v  <  f ,  x  - 


N  d  -  f  f>) 

_ D _ 

Mr  d 


[6] 


where  D  Is  tie  (interval'  icraam  of  *  d  ip  its  1  ecjr-  and  M  is  the 
value  of  f  at  its  mod°  . 

When  tils  Monte  Carto  procedure  was  used  to  generate  the  gamma  dis¬ 
tribution,  c. trrpur.atioc.  time  per  ratioi.  c.nc.r-.*  was  ie * reased  over  the  itera¬ 
tive  me* *ods  by  a  factor  of  abou*  100C.  However,  for  certain  parameter 
values,  the  gamma  distribution  is  highly  peaked;  and,  although  computation 
time  is  satisfactory,  the  method  is  not  as  efficient  as  it  could  be. 
Efficiency  is  imp or tan*  not  because  of  small  differences  in  computation 
time,  but  because  a.:,  inefficient  Monte  Carlo  program,  may  use.  a  large 
portion  of  the  power  residue  secuer. -e  to  generate  a  small  sample.  Gener¬ 
ating  large  samples  may  entail  frequent  changes  in  the  basic  uniform 
number  sequence;  tnese  charges  should,  be  avoided  whenever  possible. 

To  increase  the  efficiency  cf  the  gamma  routine,  a  modified  program 
was  written.  In  this  program,  the  x- domain  of  f  is  divided  into  sub- 
intervals  of  equal  length  For  each  subinterval,  values  of  y  are  gener¬ 
ated  only  in.  the  subset  of  the  range  of  f  which  is  supported  by  the  x 
subirterval.  This  restriction  or.  the  y  values  increases  the  efficiency, 
because  the  fa- tor  M*d  (Equation  6)  Is  replaced  by  a  sum  of  small  rec¬ 
tangular  areas  whose  total  approaches  J  f(x'ldx  as  their  number  increases. 

D 

With  even  a  small  cumber  of  intervals,  the  efficiency  is  improved  by  a 
factor  of  2  or  5.  Table  1  illustrates  the  difference  between  the  two 
methods . 


A  Monte -Carlo  method.  A  pair  of  independent  random 
numbers  determines  a  point  (x,y).  If  y  falls  in  the 
shaded  area,  x  is  discarded;  if  y  falls  beneath  the 
curve,  x  is  added  to  the  sample. 


EFFECT  OF  REF  INEMENT  ON  vcrT  CARLO  GENERATION 
OF  THE  '^AMMA  D 'JFTRIE’JTTGN* 


cairs  of  Fauactt 

Accuracy  of  Fit  to 

r 

Numbers  Needed* 

C-ama  Listributlon1 

1.0 

54.717 

*1*57 

Original 

Mel  hod 

3.0 

1.3  551 

2.19 

5.0 

13  810 

-0.35 

1.0 

6  414- 

-0.03 

Modified 

Me trod 

3.0 

6  30C 

-0.73 

5-0 

6,000 

0.34 

‘Each  row  of  the  table  represents  a  simulation  ran.  Toe  parameter  r 
(see  Equation  1 )  determines  tuft  "s  tape"  of  the  distribution.  Wien  r  is 
unity,  the  distribution  is  Eta.: ..y  skewed,  with,  a  mode  at  zero.  As  r 


Increases  toe  3  d-  (which  is  a*,  tie  Dels'  : - •  moves  to  the  right,  and 

*  >* 

A  * 

the  distribution  becomes  flatter  In  eac*  simulation,  A  was  unity,  and 
tie  same  randan?  number  sequence  was  used. 

bThe  table  entry  is  tie  number  :f  pairs  needed  to  gene T ate  a  sample  of 
5000, 

The  output  w-aa  a  100-cell  histogram ,  similar  fee  ti-e  one  srsown  In  Figure  4 
Chi-square  goodness -of • fl*  values  are  approxima'  ly  normally  distributed 
when  tie  number  of  cells  is  30  or  more.  The  table  entry  Is  toe  deviation 
In  standard  deviation  units,  of  the  computed  chi-square  value  fr'jm  Its 
theoretical  mean.  A  9^  confidence  Interval,  symmetric  about  the  theo¬ 
retical  mean,  would  have  Its  boundaries  at  ±1.96  standard  deviations. 


Composition  Methods.  Occasionally,  it  is  easier  to  convert  ore  dis¬ 
tribution  into  another  by  a  method  ether  than  integration.  This  would 
certainly  be  the  case  if  we  desired  the  gamma  distribution  only  for  inte¬ 
gral  values  of  its  parameter  r,  because  the  sum  of  r  independent  identi¬ 
cally  distributed  exponential  random  variables  has  the  required  gamma 
distribution  There  are  racy  special  methods  which  can  he  characterized 
by  tneir  reliance  on  composition  per at ions ,  such  as  addition,  multipli¬ 
cation,  or  convolution.  In  the  case  of  the  normal  distribution,  the 
central  limit  theorem,  which  guarantees  the  asymptotic  normality  of  sums 
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of  independent,  random  variables,  suggests  an  alternate  way  to  produce 
normally  distributed  numbers.  Tie  uniformly  distributed  numbers,  {*}, 

1  .  2 

0  <  x  <  1  have  mean  —  and  variance  =  I  x  -  dx.  If  12  Independent 

2  1^  Q  d 

values  of  x  are  added,  the  distribution  of  the  sum  will  rave  mean  6  and 
unit  variance  This  distribution  approximates  a  normal  list" lb ut i on  very 
closely,  even  though  only  12  random  variables  have  been  used.  Its  mai' 
drawback  Is  that,  being  bounded,  tne  fit  tc  the  normal  is  not  good  In  the 
tails.  Figure  3  shows  the  output  of  a  subroutine  which  sums  2d  uniform 
variables.  This  distribution  has  been  transformed  to  approximate  -.he  unit 

normal  bv  subtracting  12  from  each  sum  and  multiplying  the  result  by  > 

The  familiar  application  of  the  Poisson  distribution  tc  "completely 
random"  events  (e.g.,  telephone  traffic)  suggests  a  special  method  *c  con¬ 
vert  uniform  randan  numbers  to  a  Poisson  distribution.  Suppose  * i.a*  in 
a  -unit  of  time,  the  number  n  of  occurrences  of  an  event  (e.g.  telephone 
calls)  is  a  Poisson  distributed  random  variable,  l.e.. 


P(n  =  k)  =  -~r~  e  "v  ,  C7] 

where  v  >  0,  k  =  0,  1,  2,  ... 

It  follows  (Feller,  1962)  that,  If  t,,  t^,  ...  are  the  times  at  which  the 

events  occur,  tne  inter -occurrence  times,  In  =  t.,  -  t.  , 

I  =  t,  -  t_  ...  are  independently  exponentially  distributed.  Thus  ,  if 
j  3  ^ 

the  computer  generates  a  sequence  of  exponential  randan  variables, 
y.  ,  y_,  ...  ar.d  adds  them,  the  smallest  value  of  n  such  that 


2  y  >  1 

i=l 


will  be  Poisson  distributed.  This  method  can  be  simplified  by  noting 
that  by  integrating  Equation  2,  we  can  express  the  y's  in  terms  of  uniform 
numbers  > 


-In  (x, ) 
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.3 


Normal  Deviates 


Figure  3 


Output  histogram  of  the  sum  of  24  uniformly  distributed  random 
variables.  A  normal  distribution  with  zero  mean  and  unit 
variance  is  drawn  through  the  histogram,  which  represents  a 
sample  of  2000  observations. 


where  13  a  uniform  rand  our  number  and  v  is  the  parame*er  of  the  ex¬ 

ponential  distribution.  quality  8  car.  cow  be  wr  - 


or , 


v 


<  1 


and,  taking  cx£cm.ei  tie  Ls ,  V..*  smallest  values  ctf  t 


tea*. 


1-1 
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will  be  Poisson  distributed  with  mean  and  variance  v.  For  small  values 
of  v.  relatively  few  x  's  are  reeled,  and  the  routine  work?  efficiently. 

JL 

For  larger  values  of  v,  the  computation  is  slow,  tut  because  the  Poisson 
distribution  approaches  the  normal  distribution  with,  increase  In  v,  the 
Inefficient  part  of  the  routine  need  never  be  used. 


CORRELATED  A3®  NOW- IJffiE PENDE NT  RANDOM  NUMBERS 

One  requirement  of  a  simulation  program,  may  be  to  reproduce  correla¬ 
tions  that  are  observed  between  numbers  in  empirical  data.  If,  for  ex 
ample,  a  matrix  of  simulated  test  scores  is  needed,  it  is  very  often 
desirable  to  simulate  the  row  or  column  intercorrelations  as  well  as  *  he 
distribution  properties  of  the  scores-  Suppose  R  Is  the  desired  matrix 
of  inter correlations .  If  a  vector  A  of  numbers  generated  by  the  computer 
is  uncorrelated  (i.e.,  the  expected  value  of  l/n  A  A,  where  n  is  the 

i_ 

number  of  sets  of  values,  is  the  identity  matrix),  then  B  =  AR2  is  a 
transformation  of  A  which  has  the  same  distribution  properties,  but  its 
expected  intercorrelation  matrix  is 


i_  l 

B:B  =  R2"  A:  A  R2 

k-  i. 

*  R3  R2  [.Ll] 

=  R„ 
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Thus,  by  a  simple  multiplication  operat ion- -but  one  that  often  could  not 
be  done  without  a  high-speed  computer— the  random  number  output  reflects 
the  intercorrriations  of  the  empirical  da 'a.  This  method,  and  the  adequacy 
of  the  programs  using  it,  are  discussed  ir  boldt,  Wiskoff,  and  Fitch  (i960). 

The  correlation  adjustment  leaves  successive  score  vectors  statisti¬ 
cally  independent,  because  the  parameters  of  the  generating  process  are 
left  unchanged.  When  dependent  random  numbers  are  required,  the  generating 
algorithm  must  adjust  its  parameters  on  trial  (or  time  period)  t  as  a 
function  of  the  output  sequence  through  trial  t  -  1. 

Probability  mechanisms  with  time -dependent  parameters --known  as 
stochastic  processes--are  common  elements  of  simulation  models.  The 
simplest  sequential  adjustments  of  parameters  are  those  that  depend  only 
on  the  trial  numbers,  t,  and  not  on  the  magnitude  of  the  rand  cm  numbers 
themselves.  For  example,  suppose  an  exponential  random  number  generator 
is  altered  so  that  the  parameter  X  is  proportional  to  the  unfilled 
portion  of  the  desired  sample.  If  N  is  the  total  sample  size,  the  first 
value  of  X  is  X^  =  kN,  for  some  fixed  k.  On  the  second  trial,  only  N  -  1 

numbers  are  required,  and  =  k  (N  ~  l).  Continuing  in  this  way  we  have. 


Xt  =  k(N  -  t  +  l),  t  =  1,2,... N.  [12] 

Eecause  X  becomes  smaller  each  time  a  number  is  added  to  the  sample, 
successive  numbers  are  drawn  from  exponential,  distributions  (Equation  2) 
with  steadily  increasing  means  (since  the  mean  of  the  exponential  dis¬ 
tribution  is  Inversely  related  to  X).  If  we  identify  the  random  number 
output  with,  say,  response  latency,  the  generator  is  a  model  of  an 
organism  which  responds  at  a  decreasing  rate.  This  program  has  found 
useful  application  in  simulation  of  visual  search. 


CONCLUSIONS 

Tne  techniques  just  described  are  building  blocks  for  a  wide  variety 
of  randan-number  programs.  In  each  case,  the  method  has  been  programmed 
and  tested  for  efficiency  and  accuracy.  A  chi-square  testing  program  was 
written  to  evaluate  the  goodness  of  fit  of  the  generated  data  to  theo¬ 
retical  functions.  In  the  case  of  the  gamma  distribution  (which  involved 
the  most  extensive  programming  work),  typical  chi-square  values  for 
samples  of  10,000  numbers  were  significant  beyond  the  0.01  probability 
level.  This  Indicates  that  Monte  Carlo  techniques  preserve  the  basic 
accuracy  of  the  uniform  power  residue  core  program  and  can  therefore  be 
expected  to  operate  satisfactorily  in  future  applications. 
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To  evaluate  a  ra'dcc.  •  jal-r  pr«*araT  *ve  questions  trust  he  asked 
"Do  the  data  fit  t tie  de?i**i  distribute.'  and  "Are  successive  furthers 
independent?''.  All  core  ’.aerations  o;  rrccram  effieie-.cy  are  secondary 
to  these  tasic  goals,  be  aus-  :a  applicaMoe  no  aaour.t  of  computation 
speed  can  compensate  for  error,  bias  cr  ur  watted  se quest lal  dependencies, 
lie  independence  question  is  at  or  e  •  :e  most  c  ff  i  cul*  to  answer  and  the 

most,  common  source  of  *routl-  Tee  Fr  tiem  may  become  acute  when  succes¬ 

sive  random  numbers  enter  into  recursive  or  progressive  computations  *  bat 
magnify  error.  Basic allv  tee  difficulty  stems  from  systematic  variations 
in  the  uniform  core  program  Although  Tte  rov* r  residue  method  is  tie 
test  available  it  is  known  to  have  cyclic  and  coKtiensatory  props  r  *  le  s 
that  deperd  or.  t.v*  initial  number  and  multiplier  c rosea.  Toe  efi>cf  of 
these  regularities  In  a  sequence  that  ought  to  be  irregular  Is  no*  always 
easy  to  detect.  Figure  4  shows  one  instance  ween  toe  malfunction  was 

discovered  immediately,  because  the  ou*put  histogram  was  <'om.poped  of 

relatively  small,  class  interval?  *  n-  :.;*--vai  sue  had  tee:  d'ubled 

the  fault  could  Lave  escaped  notice.  Jr  principle,  the  mat hena~cr o  of 
number  theory  could  be  used  to  discover  which  choices  of  initial  r  •  ai 
t.ions  produce  the  "best"  random  sequence,  but  the  work  would  be  sc  cm- 
pi  ex  ard  tedious  that  nc  one  is  likely  to  at  temp*  it.  The  only  recourse 
is  to  trial  and  error.  Attempts  to  reduce  sequert lal  effects  by  composi¬ 
tion  of  independent  generators,  rave  been,  successful  ,  v  Tar:'  at i 
Marsaglia  196b ) . 

Occasionally,  random  numbers  are  needed  ir.  large  quantities  rather 
*ran  one  >y  one  and  the  principal  requirexint  is  that  ‘•-ev  closely 
approximate  a  distribution  function.  Some  i *.*•■? grat  ion  problems  have  *  his 
requirement.  The  independence  of  successive  numbers  In  the  sequence  can 
be  sacrificed  if  true  to*al  sample  distribution  can  be  made  precis-.  If 
the  desired  distribution  is  uniform  cr  is  based  or.  a:,  underlying  u* if on 
sample.,  the  output  of  the  random  number  generator  car  be  improved  it  this 
respect.  For  example,  if  a  sample  of'  1000  numbers  is  desired  in  a  uniform 
distribution  ever  the  uni*,  interval,  it.  car  be  achieved  by  selecting  each 
number  in  the  list  0.001,  0.002, .. .1.000  exactly  once.  A  program,  called 
RANFEF,  las  been  written  to  make  such  selections  ir.  a  random  order. 
Briefly,  this  program  uses  the  power  residue  subroutine  to  randomly  select 
the  numbers,  one  at  a  time;  each  selection  is  then  removed  from  the  list 
of  available  numbers  and  the  final  output  is  a  permutation  of  the  list . 
Because  the  numbers  are  selected  without  replacement,  the  selection  ect 
abilities  change  as  the  list-  diminishes.  Consequently,  the  output  no 
longer  random  it  is  a  haphazardly  constructed,  hut.  perfectly  uniform 
sample . 

The  divers*  applications  of  random  number  generators  make  it  likely 
that  rew  techniques  will  be  found  and  put  to  use  in  future  work..  To  da**, 
Vere  have  been  more  suggested  appl  1  cations  for  these  nrograms  than  the 
organization  has  resources  t  o  our  sue.  Toe  present  memorandum  thus  reports 
only  the  beginning  of  a  growing  research  field.  Whether  the  rew  develop¬ 
ments  wll.1  simply  b<-  more  efficient  ’•efinements  cf  basic  technique?  or 
fundamentally  different  approaches  u*.  qu-*s*ton  that.  will  he  answered 
only  with  eontiauid  iav«®tigaMar... 


Figure  4.  Output  histogram  from  the  gamma  distribution  generator. 

Each  vertical  line  represents  the  sample  proportion  in  the  ith 

class  Interval.  The  parameters  for  this  run  were  \«1,  r-3>  and 
the  sample  size  was  5000.  The  histogram  Bhows  the  effect  of  a 
pronounced  regularity  in  the  power  residue  computations  for  the 
basic  uniformly  distributed  random  numbers:  the  sample  proportions 
should  conform  to  a  function  that  lies  between  the  tall  and  short 
segments. 
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